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Abstract This paper reports an application of gas-kinetic BGK scheme to the computation 
of turbulent compressible convection in the stellar interior. After incorporating the Sub- 
<~i , grid Scale (SGS) turbulence model into the BGK scheme, we tested the effects of numer- 

Qh' ical parameters on the quantitative relationships among the thermodynamic variables, their 

Q I fluctuations and correlations in a very deep, initially gravity-stratified stellar atmosphere. 

;h ' Comparison indicates that the thermal properties and dynamic properties are dominated by 

lyj I different aspects of numerical models separately. An adjustable Deardorff constant in the 

Ci . SGS model c^ = 0.25 and an amplitude of artificial viscosity in the gas-kinetic BGK scheme 

6*2 = are appropriate for current study. We also calculated the density-weighted auto- and 
cross-correlation functions in Xiong's ( 119771 ) turbulent stellar convection theories based on 
^ ' which the gradient type of models of the non-local transport and the anisotropy of the turbu- 

^O ! lence are preliminarily studied. No universal relations or constant parameters were found for 

these models. 
I> 
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1 INTRODUCTION 

Turbulent convection has tight relations to the unsolved problems in the theory of stellar structure and evo- 
lution, especially, for the massive stars (Deng et al. 1 1 996al 1 1 996bl and references therein). These problems 
cannot be settled with pure analytical ways. Along with the development of computing science, numerical 
simulations become a powerful tool to investigate the hydrodynamic properties of astrophysical flows. It is 
widely used in the study of formation of cluster, accretion disk and evolution of galaxies. The convection in 
stellar interior have also been studied by many authors with numerical experiments. Due to the difficulties 
in this problem, the progress made in this field is limited. However, it is generally believed that numeri- 
cal testing of some analytical models and local high-resolution simulations can make a lot of sense to our 
understanding of stellar convection. So far, the numerical hydrodynamic scheme applied to the stellar con- 
vection are Lax-Wendroff scheme (Graham |1975l etc.), alternating direction implicit method on staggered 
mesh (ADISM) (Chan et al. TWT; 1986' etc.), pseudo-spectrum scheme (Hossain & MuU an [19901 [19911 
etc.), piecewise-parabolic method (PPM) (Porter & Woodward 1994; 2000, etc.), upwind scheme and so 
on. At the present time, the most suitable numerical scheme for the turbulent flow is the spectrum method, 
but it cannot handle the discontinuity in the motion of fluids. 

Gas-kinetic BGK scheme is a recently matured method for computational fluid dynamics and attached 
much attention in the practical problems (Xu 1200 11 1. It is accurate and robust for computing supersonic 
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unsteady flows. However its application to astrophysical flows is not popular attributed to its complex. 
Theoretical analysis shows that near the surface of stellar envelope, the motion of fluids becomes supersonic 
which cannot be self-consistently treated by traditional mixing-length theory (MLT) (Deng & Xiong|2001j. 
Our original aim is to simulate the supersonic turbulent convection in the outer region of yellow giants by 
gas-kinetic BGK scheme which would involve many efforts in different directions. We have already extend 
the BGK scheme to include the gravitational acceleration (Tian et al. '2007'). Before using it to compute 
the supersonic stellar convection, the turbulence model, radiation transfer model and realistic input physics 
must be correctly implemented. 

Restricted by the capacity of digital computer, we cannot afford the very high resolution numerical 
experiments for the stellar type of turbulent convection. Large eddy simulation (LES) which calculates the 
large eddies explicitly while mimics the sub-grid eddies by models may be the most feasible way in current 
stage. In current paper we implement the SGS turbulence model (Smagorinsky 1963; Deardorff 1971) into 
the BGK scheme and validate the three-dimensional BGK code by calculating the turbulent compressible 
convection in a deep stellar atmosphere. For very high Reynolds number, the behaviors of the turbulent 
flows are greatly affected by the numerical and physical dissipation in the scheme. A investigation of these 
effects is very necessary before the code is applied to the practice. By varying the Deardorff number in 
the SGS model and artificial viscosity parameter introduced to capture the shock, we constructed three 
models which are similar to those studied by Chan & Sofia (.1989 hereafter CS89; |1996l hereafter CS96). 
The empirical relations derived by them were re-examined. A study of density-weighted auto- and cross- 
correlation function, anisotropy of turbulence and diffusive type of models of non-local turbulent transports 
in the turbulent stellar convection theory of Xiong (1 19771 1 was conducted too. 

In the next section, we give a description of gas-kinetic BGK scheme and mainly focus on the incor- 
poration of SGS model. The computed physical models are formulated in Sect. [3] The numerical results 
are shown in Sect. |4] where the discussions are also presented. The conclusions are summarized in the last 
section. 

2 GAS-KINETIC BGK SCHEME 

General numerical method for hydrodynamic problems is to directly discretize the Navier-Stokes equations, 

dp/dt = -V ■ pv, (1) 

dpv/dt = -V ■ pvv -Vp + V ■Y. + pg, (2) 

dE/dt = -V ■ [{E + p)v - V ■ -S + Fd] + pv ■ g, (3) 

where p is the density, v is the velocity, p is the pressure, g is the gravitational acceleration and E is the 
summation of internal energy and kinetic energy. 

S = 2p(T + <,{\7 -v)! 

is the viscous stress tensor, where cr is the strain rate tensor, I is the identity tensor, p and <r are the 
dynamical and bulk viscosity coefficient, respectively. Fd is the diffusive type of energy flux. Differently, 
the BGK scheme works on the BGK equation, 

-± + c-\7f + g-\7J=^ L, (4) 

Ot T 

which is an approximation of Boltzmann equation (Bhatnagar et alll9541). In above expression, f{x, c, t) 
is the gas distribution function in the phase space, c is the particle velocity, t is the collision time, and 
Vc = ((9/9ci, d/dc2, d/dc^). The right-hand side of equation (|4|) is the so-called relaxation model, which 
is a simplification of the complicated collision term in the Botlzmann equation (Vincent et al. 1965). It 
physically means that the initially non-equilibrium distribution / will approach the equilibrium state f^'^ 
after the particles collide once. A larger r corresponds to a further state / from Z"^"?, i.e., the stronger non- 
equilibrium transport effects, such as viscosity and conduction. In our study, the equilibrium state, f'^'^, 
in equation dUi is taken to be MaxwelUan distribution. It can be proved mathematically that the solutions 
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of Navier-Stokes equations ([TJ-O are automatically obtained through solving the BGK equation with the 
following definitions of dissipative coefficients: 

2 A/" A/" + 5 fc 

fi = Tp, <;=--— —Tp, K= — rp, (5) 

3 A/ + 3 2 m 

where k is the thermal conductivity, J\f is the internal degree of freedom for particles, m is the molecule 
mass and k is the Boltzmann constant. 

The basic idea of BGK method is to find a local approximation to the non-linear equation (|4]i. Then 
evaluate the macro quantities (e.g., fluxes) using the micro distribution function /. Finally, the cell average 
values are updated according to the conservation laws (finite volume method). In our first attempt, the 
BGK method (Xu 1200 il l has been extended to include the external force. A detailed description of the three 
dimensional multidimensional gas-kinetic BGK scheme for the Navier-Stokes equations under gravitational 
fields was given by Tian et al. (2007) . Here, we only outline the new implements, i.e., the incorporation of 
turbulence model. 

In the study of convection, the combined effects of viscosity, heat conduction and temperature differ- 
ence on the instability of flows are measured by Rayleigh number: Ra — galSTd^ /{vn), where a is the 
thermal expansion coefficient and v is the kinematic viscosity. For the polytropic gas defined in Sect. |3] the 
Rayleigh number can be written as: 



Ra = 



PrRTtZ^d^pi 



^? 



1 - (7 - l)n 

7 



(" + !), (6) 



where R is the gas constant, d is the depth of computational domain and 7 is the ratio of specific heat, the 
meaning of other symbols can be found in Sect. [3] In the gas-kinetic BGK scheme, the viscosity is controlled 
by collisions of particles. We can relate the Rayleigh number to the collision time by the following way. 
Suppose during each time-step Ai, the particles collide (3 times. Then we have 

A^ = /5r=^, (7) 

where 5 is Courant number. Ax is the spatial resolution, v is speed of fluids and Cg is the sound speed. From 
equation (|5]l, ^ and (|7|l, we get 

Ra^ ^ ^[l-(7-lH(7i + l)/32, (8) 

where Ma = vjcg is the Mach number, N is the vertical grids size and Pr is Prandtl number. For current 
study, we have N ^ 50, Ma ^ 1, Z = 15, (5 = 0.3, 7 = 5/3 and n = 0.999/(7 - !)■ In efficient turbulent 
convection, the energy transfer by heat conduction is negligible which means the Prandtl number is very 
large. In all of our simulations, the Pr is set to be 10^. Therefore, approximately we have max (Ra) ~ 
6 X 10^/32. Similarly, we have Reynolds number Re = NMail + Ma)(il5 - 3 x W[3. While in flie 
typical stellar convection. Re ^ 10^". Hence, /? ^ 10® is needed when we do a direct numerical simulation 
(DNS) of stellar convective flows where Ra ^ 10^"*. These values can be reduced by increasing grids 
number which is very expensive for the present generation hardware. At the same time, a very small r 
would introduce large computational error. An alternate way is the LES which simulate the large eddies 
directly and approximate the small eddies with models. There are a lot of approaches to perform the LES, 
the simplest way may be the SGS model (Smagorinskv l 19631 Deardorff |I971| ). 

In the BGK scheme, the viscosity is introduced through the collisions of particles. The natural way of 
implementing SGS model is to modify the collision time. Chen et al. ( I2003I I included the renormalization 
group k — e large eddy model into the BGK equation, where k is the turbulent kinetic energy and e is the 
turbulent dissipation. In the fc ~ e model, tow additional equations are needed to be solved. For sake of 
simplicity, we consider the Smagorinsky ( |I963l l model. In current study, the collision time is defined as 

Ttot = T + ij2-. ; rAtH , (9) 

\Pl+Pr\ P 
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where the first term of right-hand side represents the molecule viscosity and the second term is introduced 
to increase the numerical dissipation when there is a jump in pressure around the control volume boundary. 
C2 is an adjustable constant, pi and pr are the reconstructed pressures at the left and right side of a cell 
interface (see sketch in Fig. 1). In the strong supersonic region, additional dissipation caused by this term is 
essential to stabilize the computing. 



i+1 



p,pv,E,p) 



p,pv,E,p) 



Fig. 1 Reconstruction of conservative variables and pressure where discontinuity is introduced 
at the cell interface. The dashed lines represent the boundaries of the of the control volumes 
numbered by i and j; the vertical solid line is the interface where the flux-splitting is performed. 
The discontinuity of the values and their slopes are depicted by oblique solid lines. 



The last term in right-hand side of equation (|9]l is implemented to account for the SGS viscosity. 



f-^sgs 



p{c^Af{2(T:(T[ 



,1/2 



(10) 



where Cn 



u^ is an adjustable constant, usually has the value: 0. 1 
resolution, colon stands for the contract of tensor and cr = diV 



0.2, the filter width A is taken to be the local 
J +djVi. Above model is called Smagorinsky 
model or sometimes Smagoringsky-Lilly model. In our calculations, the eddy viscosity is computed in the 
control volume by staggered mesh strategy and then interpolated to the cell interface. In the BGK scheme, 
there is intrinsic diffusion caused by particle collisions. In current study, we are only interested in the 
turbulent properties. So, the molecule Prandtl number Pr is set very large. And the following diffusive flux 
(CS96) is implemented explicitly into the BGK scheme. 



-Ct'^T - CsVS, 



(11) 



where S = Cp{\nT — Vq Inp) is the specific entropy, Cp is the specific heat at constant pressure and Va 
is the adiabatic gradient. In the stable layer, Ct is set to make the diffusion carry out the input energy flux. 
In the convection zone, Ct is very close to zero. Cg = Hsgs/Pfsgs represents the turbulent diffusion and 
is set to zero in the stale region, where Prggs is the effective Prandtl number of SGS turbulence and taken 
to be 1/3. 

3 PHYSICAL MODELS 



Our physical problem is very similar to those studied by CS89 and CS96. An ideal gas (p — pRT) in a 
rectangular box is considered with gravity in the vertical direction. The side boundaries are periodic and the 
top and bottom boundaries are impenetrable and stress free. In order to avoid boundary effects, a very thin 
stable layer is placed below the upper boundary and the diffusive flux is gradually enhanced near the lower 
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boundary to make it carry out total flux at the lower boundary. A constant energy flux Ff, = 0.25 is fed at 
the bottom. At the top, the entropy is fixed. The system is initially static: 

(1 + Z{d - z)/d)Tt, (12) 

{T/Ttrpt, (13) 

(T/T,)"+Vt, (14) 

{n+l)ptZ/{ptd) (15) 

where Z — {Tb — Tt)/Tt is the normalized parameter, < 2; < 1 and n is the polytropic gas index. 
The gravitational acceleration g comes from the hydrostatic equilibrium dp/dz = —pg- The subscripts t 
and h denote top and bottom values respectively. Above solutions to Navier-Stokes equations are not stable 
against small perturbations. In all of our calculations, the velocity field is slightly perturbed initially. After 
long-time thermodynamic relaxation, the system will reach a statistical steady state. We defined a series 
of runs to test the effects of numerical parameters. The numerical effects of a variety of parameters were 
studied by Chan & Sofia (119861 hereafter CS86) and CS89 in details. The effects of changing turbulent 
Prandtl number was tested by Singh & Chan (ll993l l. Here, we just focus on the important parameter c^ in 
the SGS model and new parameter C2 appearing in the BGK scheme. The details are given in the second 
line of Table [T] 

All the cases we computed use a 29 x 29 x 45 mesh. The vertical grid decreases smoothly with height 
(about 6 grids per PSH) and the horizontal grid is uniform. The aspect ratio (width/depth) of the box is 1.5. 

4 RESULTS 
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Fig. 2 Distribution of temporally and horizontally averaged energy fluxes with depth. Dash dot 
line: diffusive flux; dotted line: kinetic energy flux; diamonds: enthalpy flux; solid line: total 
energy flux. 



In this part, we show the results from the numerical simulations. All the runs were evolved 2000000 
numerical time-steps, corresponding to a dimensionless time around 858, before the statistical analysis is 
performed. The statistical steady state is indicated by the balance of the input energy flux from the bottom 
and the outgoing energy flux through the top. In our calculations, the spatial variation of averaged total 
energy flux from 0.25 is within 0.1% (see solid line in Fig. 2). Another criterion is the averaged vertical 
mass flux which is less then 10^^ everywhere in all cases. Hence, the system will not undergo substantial 
adjustment any more. The statistical average covers 500000 numerical time-steps. 

Except for the instant velocity fields, all the other quantities investigated here are the mean values. For 
an arbitrary quantity q, (q) represents its combined horizontal and temporal mean, q' denotes the deviation 
from (q), q" stands for the root mean square (rms) fluctuation from the (q). 
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In some figures, the integral pressure scale height (PSH) is shown by the vertical lines. For example, in 
Fig. 2, the vertical solid line denotes the location of stable-unstable interface near the upper boundary. The 
second dashed line at the left side of the solid line is 2 PSHs away from the upper stable-unstable interface. 
Although the numerical parameters are different for case A to case C, their relaxed thermal structure are 
nearly equivalent and the discrepancy between the integral PSH locations is very small. So, we can plot the 
results from all runs in the same figures. 

4.1 Velocity Fields 





Fig. 3 Instant velocity fields at time 858. Left: projection in vertical plane at y = 0.45; right: 
projection in horizontal plane at z = 0.6. 



In CS86, the three-dimensional turbulent flow structure was well depicted by the pseudo stream lines. 
For sake of clarity. Fig. 3 shows the velocity fields projected in x-z plane and x-y plane. From the left 
panel of Fig. 3, it is evident to see that the high speed motions exist in the top region and are associated 
with the downward streams. In our calculations, the maximum Mach number Ma occasionally exceeds 



one. In CS96, the SGS viscosity was enhanced by a factor [1 + 0.5(w^ 



v^)/c^] to suppress the 



shocks occurring in the top region which would easily trigger the instability of the numerical computation, 
especially, during the early thermal relaxation. Based on flux-splitting (see Fig. 1), such implement can be 
avoided in gas-kinetic BGK scheme. If the nonlinear van Leer limiter is replaced by central interpolation, 
the supersonic motion cannot be handled correctly. Numerical tests also show that a Ra less than 10^ will 
smooth the turbulence, the circulations become laminar The solutions are almost unchanged when Ra great 
than 10^^. In current study, we adopt Ra ~ 10^^. The networks of downward streams can be seen the right 
panel of Fig. 3. 



4.2 Approximate Relations 



Our BGK code has been extensively tested for laminar flows (Tian et al. I2007I I. In order to validate the 
incorporation of SGS model, we re-estimate quantitatively some approximate relations among the thermo- 
dynamic variables and their fluctuations. The results are given in Table[T]where the results from CS89 are 
also listed for comparison. Our computational models are partially similar to CS89 and partially to CS96. 
Models of CS89 undergo substantially adjustment near the boundaries and we cannot afford the high- 
resolution desired in CS96. During the data analysis, the correlation function of quantity q and p is defined 
as C{q,p] = (qp) / {{q'^y ^'^ (p'^y ^^) ■ The standard deviations of these approximations (ai) are given in the 
brackets. The deviation of i?i fromCS89 is calculated by (72 = \{{Ri A+Ri B+Ric) / ^— Ric S89) / R^c SS9\- 
We only concentrate on the middle region in the convection zone, namely, 1 PSH from the bottom and 2 
PSHs from the upper stable-unstable interface. The investigated layer expands about 3 PSHs. The approx- 



Table 1 Quantitative Estimates of Some Approximate Relations. 



Relation 



Parameter 



CS89 



Deviation ID Category 



Cm 



{P"/{P))/{T"/{T)) 
{p"/{p))/{T"/{T)) 
S"/{CpT"/{T)) 

{p"Kp))IW'V{t)) 
v"KpW 

CiT',S') 
C{p',S') 
C{p',T') 
C{p',T') 

C{v.,T') 

C{v.,S') 

C{v.,p') 

{v.p'}/{v.){p) 

{v,p)/{v,){p) 

{v.T')/{v.){T) 

{v.S')/Cp{v.) 

MKv"l){T) 

Fep/Cp{p){v,) 

F,p/Cp{p)v"l 
T"/{T) = aAV + b 



v"/{T) = aAV + b 



0.2 
1 



0.25 
1 



0.2 




0.2 



(%) 



AV = a[F6/(0.8Cp(p){T)i/2)]2/3 + ^ 



0.69(0.12) 

0.57(0.10) 

0.83 ( 0.01 ) 

0.57 ( 0.04 ) 

0.90 ( 0.01 ) 

0.48 ( 0.01 ) 

0.88(0.14) 

1.52(0.15) 

0.98 ( 0.00 ) 

-0.93 ( 0.01 ) 

-0.83 ( 0.02 ) 

0.59 ( 0.04 ) 

0.78 ( 0.02 ) 

0.76 ( 0.02 ) 

-0.64 ( 0.01 ) 

-1.00(0.02) 

1.49 ( 0.06 ) 

1.46 ( 0.06 ) 

1.28(0.05) 

0.81 ( 0.08 ) 

1.49 ( 0.06 ) 

1.20(0.15) 

a= 1.92 

b = 0.0029 

( 0.0004* ) 

a = 1.01 

b = 0.0027 

( 0.0004* ) 

a = 0.71 

6= -0.0013 

( 0.0002* ) 



0.66(0.11) 

0.58(0.10) 

0.83 ( 0.02 ) 

0.56 ( 0.04 ) 

0.90 ( 0.01 ) 

0.47 ( 0.01 ) 

0.85(0.13) 

1.51 (0.13) 

0.98 ( 0.02 ) 

-0.92 ( 0.01 ) 

-0.83 ( 0.02 ) 

0.59 ( 0.04 ) 

0.79 ( 0.02 ) 

0.78 ( 0.03 ) 

-0.67 ( 0.01 ) 

-0.99 ( 0.02 ) 

1.44 ( 0.05 ) 

1.41 ( 0.06 ) 

1.24 ( 0.04 ) 

0.85 ( 0.07 ) 

1.44 ( 0.05 ) 

1.23(0.13) 

a = 1.88 

b = 0.0025 

(0.0001* ) 

a = 1.03 

b = 0.0023 

( 0.0002* ) 

a = 0.71 

b= -0.0011 

( 0.0002* ) 



0.62(0.11) 

0.64(0.11) 

0.83(0.01) 

0.57 ( 0.04 ) 

0.90 ( 0.01 ) 

0.48 ( 0.00 ) 

0.87(0.14) 

1.51 (0.15) 

0.98 ( 0.03 ) 

-0.92 ( 0.01 ) 

-0.82 ( 0.03 ) 

0.59 ( 0.03 ) 

0.77 ( 0.02 ) 

0.75 ( 0.02 ) 

-0.63 ( 0.01 ) 

-1.00(0.02) 

1.50(0.06) 

1.46 ( 0.06 ) 

1.28 ( 0.05 ) 

0.80 ( 0.09 ) 

1.50(0.06) 

1.20(0.15) 

a = 1.95 

b = 0.0030 

( 0.0002* ) 

a = 1.03 

b = 0.0028 

( 0.0004* ) 

a = 0.70 

b = -0.0013 

(0.0001* ) 



0.61 ( 0.05 ) 

0.61 ( 0.05 ) 

0.89 ( 0.04 ) 

0.57 ( 0.07 ) 

0.94 ( 0.03 ) 

0.26 ( 0.01 ) 

0.51 (0.03) 

0.90(0.10) 

0.99 ( 0.01 ) 

-0.89 ( 0.06 ) 

-0.82 ( 0.05 ) 

0.49 ( 0.05 ) 

0.81 (0.03) 

0.81 (0.03) 

-0.74 ( 0.03 ) 

-1.00(---) 

1.24(0.08) 

1.26(0.08) 

1.20(0.08) 

0.58 ( 0.07 ) 

1.25 ( 0.08 ) 

0.72 ( 0.06 ) 

a= 1.05 

b = 0.0027 

( 0.0008* ) 

a= 1.17 

b = 0.0032 

( 0.0008* ) 

a = 0.9 

b = -0.002 

( 0.0008* ) 
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Fig. 4 Height distributions of fluctuations and relative fluctuations for case C. Left: p" / (p) (solid 
line), T"/{T) (dotted line), p"/{p) (dashed Une); right: w^' (dashed Une), v'^ (dotted Une), v'^ 
(solid line). 



imate relations and correlations in Table [T] are classified roughly into four categories which are indicated 
by Roman numerals in the last column according to the goodness of fit and discrepancy from CS89. These 
categories are: Category I: cti < 0.6 and fT2 < 10%; Category II: ai < 0.6 and (72 > 10%; Category III: 
cTi > 0.6 and era < 10%; Category IV: cti > 0.6 and (T2 > 10%. 

Category I contains the best relations where the thermal variables are mainly involved, i.e. p, p, T 
and their fluctuations. These relations should be weekly affected by numerical scheme and details of the 
computational models. In the work of Kim et al. (1995) where a realistic equation of state (EOS) was used, 
these relations are different from current studies. Therefore, Category I may be dominantly determined by 
EOS. 

Most relations of Category II are functions of pressure p , velocity v and their fluctuations. The ampli- 
tude of the fluctuations of the components of velocity in our computation is lower than those in CS89 (Fig. 1 
therein) and CS96 (Fig. 2 therein) while the amplitude of the relative fluctuations of thermal quantities are 
around two times larger than those in CS89 (see Fig. 2 therein) and in Singh's work ( 119931 1. Hence when the 
fluctuations of thermal variables are expressed in terms of v" , w", etc., the coefficients are nearly doubled, 
e.g., Rb and i?22. This can also be seen in Fig. 3a in CS96. 

Ratios of f "/f " and v''/v'l do not deviate from CS89 very much, but the fitting is not accurate. This 
may be caused by the effects of the boundary and transition layers. It is obvious from Fig. 4 that the height 
distribution of w" is not flat as CS89. There is a small hump below the upper transition layer (one PHS 
from the unstable-stable interface) which can also be found in CS96. We believe that the unstable-stable 
transition is responsible for this. The large difference in Category IV should be caused by the high order 
powers of velocity fluctuations. 

i?22 ~ i?24 are commonly approximated in MLT. In current study, these relations are fitted by linear 
approximation very well. However, the slopes are different from the CS86 and CS96. The situation of i?22 
can be explained as above. For i?23, the reasons may lie in the fact that in the upper convective region, the 
amplitude of w" from CS96 is larger than ours and therefor we need a smaller i?23. CS96 got a smaller i?24 
(0.78) than CS89. In current study, it is even smaller which implies a smaller super-adiabatic gradient AV. 

The cause of most of the discrepancies between current study and CS89 may be the larger T" and 
smaller w" . In the traditional theory of stellar convection, the enthalpy flux is proportional to the fluctuations 
of temperature and velocity, i.e., Fg ^ Cpv"T" and the kinetic flux is totally neglected. When the kinetic 
flux is comparable to the total flux, this kind of proportional relation is not exactly held any more. Based on 
the following facts: 



Fb = 0.0625, max T" 
Fh = 0.1250, max T" 



0.022, max < 
0.023, max v" 



0.2 (Singh [T993l Fig. 1 and Fig. 5); 
0.22 (CS89, Fig. land Fig. 2); 



3. Fb = 0.2500, max T" 

4. Fb = 0.2500, max T" 

5. Fb = 0.6551, max T" 
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' 0.065, max<' - 0.35 (CS96, Fig. 2); 

' 0.078, maxw^' - 0.23 (curr ent. Fi g. 4); 

' 0.129, max<' - 0.27 (Kim [T995l Fig. 7 and Fig. 8), 



we can conclude that the T" and u" are proportional to Fb in a nonlinear way. Note that current study use 
a different numerical scheme and Kim adopted a realistic EOS. It is beyond the scope of current study to 
perform a quantitative analysis of such relations. 

Above comparison shows that the dynamical properties and thermal properties for stellar type of con- 
vection may be affected by different aspects of the numerical models separately. Thermal structure is mainly 
determined by physical parameters while dynamic motions can easily affected by numerical parameters. 




0.4 0.6 

HEIGHT FROM THE BOTTOM 



Fig. 5 Height distributions of natural logarithm of x, V and Z. Solid line: case A; dotted line: 
case B; dashed line: case C. 



4.3 Anisotropic Turbulence 

Xiong's non-local time-dependent stellar convection theory is based on Reynolds stress method. It is dy- 
namic theory of auto- and cross-correlation functions of turbulent velocity and temperature fluctuations. 
These fluctuations are defined as the derivation from the density weighted average, i.e.. 



(pvt) 



f' ^T 



(PT) 



' ' ip)' ^ ^ ip) 
The start point of Xiong's theory is a set of partial differential equations of 

X^ = (z«>'^)/3, Z=(f'2)/(f)2, V={f'w';)/{f), 



(16) 



(17) 



where w'^ — pu[/ {p), (T) = {pT)/{p) and the summation convention is adopted. The numerical results 
of X, Z and V are given in Fig. 5. In the closure models of Xiong's theory, three adjustable parameters, 
i.e., ci, C2 and C3 are introduced. C3 is used to describe the anisotropic turbulent motions. In Deng's (I2006I I 

(3 + C3)/2c3 in the fully unstable zone. In 



work, C3 was related to turbulent velocity by w' ^/{w'^ 



the upper overshooting region, they proposed that w' ^/{w' 



Wy) 

2 



/2\ 



0.5 and is independent of C3. In the 



' w'y) < 0.5 and decrease as C3 decrease. There is no enough room for 



lower overshooting zone, w' ^/{w' ^ 

overshooting in our models. The ratio ofw'^/ {w' ^+w' ) from numerical simulations is given in Fig. 6 from 
which we can see this ratio is slightly dependent on c^ and C2 . In the upper efficient-inefficient convection 
interface (about 1 PSH from unstable-stable interface), this value approximately equals to 0.5. It takes its 
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maximum at the location where the turbulent convection starts to become inefficient near the bottom. Its 
maximum is affected evidently by c^ and C2. Current study cannot give a definite solution for anisotropic 
turbulence. Here, we present preliminary suggestions. Suppose that w' ^/{w' ^ + w') = (3 + C3)/2c3 is 
held in current models, we can conclude that C3 is infinitely large at the upper efficient-inefficient interface 
and decreases as the distance from the boundaries increases. The minimum of C3 is about 0.75. 




0.0 0,2 0.4 0,6 0.8 
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Fig. 6 Height distribution of w' ^/{w' ^ + w' ). Solid line: case A; dotted line: case B; dashed 
line: case C. 



4.4 Non-local Transport Models 

Generally, a Reynolds stress method suffers the so-called closure problem. CS96 numerically studied the 
popular closures and found they were poor. Some third order moments representing the non-local transport 
effects were approximated by gradient models in Xiong's work ( 119891 ), i.e.. 



NLTl = 
NLT2 = 
NLTi = 



{uiw'^f')/{f) = -xk^,{{w'^f')/{f)), 



(18) 
(19) 
(20) 



with Zi ~ /3 ~ /s = A, where A is the Lagrangian integral length scale of turbulence. The non-local 
transports and coefficients from numerical simulations for all cases are shown in Fig. 7 . From panel (a) of 



Fig. 7, we can see that the | {z 



is about one order larger than |(w'j,w"T')/(T}| and |(u'j.r'2}/(r)2| is 



three orders less than \{u'j.w''^T') / {T)\. Hence, the non-local transports are dominated by turbulent kinetic 
energy. Panel (b)~(d) of Fig. 7 show that the /i, /3 and ^5 gradually increase with the distance from the top 
boundary and change this trend near the interface where the turbulent convection become inefficient. The 
variation ofli is slow in the lower half convection zone around the value of 1 .8. But it is affected by c^ and 
C2 obviously. Z5 is about ten times larger than l^. Both of l^ and l^ vary rapidly which suggests that there 
may not exist the universal constant for these closure models. 



4,5 Effects of Numerical Parameters: c^ and C2 

In current LES, the local effective grid Reynolds number is enlarged by SGS model whose amplitude is 
controlled by Deardorff constant c^. An inadequately small c^ would cause the building-up of the kinetic 
energy at the two-grid level and make the computation crashed. If c^ is too large, the turbulent motions 
will be damped down. The proper value of c^ is resolution- and method-dependent since the numerical 
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Fig. 7 (a) Non-local turbulent transports (in common logarithm), (b)^(d) The approximate co- 
efficients. Solid line; case A; dotted line; case B; dashed line; case C. 



dissipation also plays an important role in the behavior of turbulence. Deardorff ( |1971| l suggested a value 
of 0.2 for c^ which was used in the CS89. In CS96, this value was increased to 0.25. 

In the flux-splitting method, the discontinuity is introduced at the cell interface (see Fig. 1) by limiters. 
Additional dissipation (the second term in the left hand side of equation (|9|l) is employed to handle the 
strong shock waves near the sharp jumps of pressure. The typical value of C2 is 1 . In the smooth region 
this kind of discontinuity should be very small. However, our simulations have lower resolution and the 
computing zone extends about seven PSHs. So it is necessary to check if these values are adequate for 
studying the stellar type of convection. 

Relations R2 ^ R4 and R9 ^ i?ll in Table [T] show that the effects of changing these parameters on the 
thermal structure are really slight. They both affect the eddy properties with very small amplitude which can 
be seen from Fig. 8 where the auto-correlations of the vertical velocity are shown. Their profiles are nearly 
symmetrical except in the upper stable and transition zone. The half width at half maximum (HWHM) of 
these profiles should be sensitive on the viscosity (CS86). RO and Rl in Table[T]are the f "/f " and v''/v'l, 
respectively. They should be nearly equal to each other for isotropic turbulent flows. In our study, for case 
A (cp — 0.2, C2 = 1), the discrepancy between them is clear which becomes slight for case B (c^ — 0.25, 
C2 ~ 1) and further slight for case C (c^ = 0.2, C2 — 0). This kind of anisotropy comes from the initially 
perturbation and can also bee found in Fig. 7 and Fig. 11 in CS86. So it seems that c^ = 0.25 is more suitable 
for current study. In the deep stellar convection, the shock waves are mild. Hence, the flux-splitting may 
enough to handle it as in our tests. We suggest that it is better to make C2 as small as possible because it 
would enhance the anisotropy which can be diminished by larger c^. 
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Fig. 8 Comparison of auto-correlations of the vertical velocity. Solid line: case A; dotted line: 
case B; dashed line: case C. -pi is the pressure at the unstable-stable interface. 



5 CONCLUSIONS 

In this paper, we present a preliminary application of gas-kinetic BGK scheme to the simulation of turbulent 
convection in stellar atmosphere. The approximate relations among thermodynamic variables, their fluctu- 
ations and correlations were examined. The anisotropy and diffusive models of non-local transport were 
investigated too. The effects of varying numerical parameters were also tested. The main conclusions are 
summarized as follows: 

1 . The behavior of the thermal variables and dynamic variables are affected by different aspects of the 
models and numerical scheme. For example, the fluctuations of density and pressure are dominantly 
determined by physical models while the fluctuations of velocity are sensitively dependent on numerical 
parameters, e.g., c^ and C2. 

2. There is no constant ratio of w"/(w" + w") for anisotropic turbulence in current models. We suggest 
that C3 take an infinite value at the boundary and approach its minimum (0.75) in the deep convective 
region. 

3. The diffusive models for non-local transport are not applicable since the coefficients for different quan- 
tities are dramatically different. The best situation is the turbulent transport of turbulent kinetic energy 
where a roughly flat region exists for l^- 

4. For current resolutions, c^ = 0.25 is better than c^ = 0.2. And C2 should be set as small as possible in 
any cases. A flux-splitting technique is needed to stabilized the shock waves near the top. A Rayleigh 
number less than 10^ may smear the turbulent motions. 

Our simulations may suffer from the lower resolutions, aspect ratio and the number of testing cases. 
But it is enough for the purpose of validation and get some preliminary results. A further study will be 
performed in the future. 
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